Phase diagram of a hard-sphere system in a quenched random potential: a numerical 

study 
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We report numerical results for the phase diagram in the 
density-disorder plane of a hard sphere system in the presence 
of quenched, random, pinning disorder. Local minima of a 
discretized version of the Ramakrishnan-Yussouff free energy 
functional are located numerically and their relative stability 
is studied as a function of the density and the strength of dis- 
order. Regions in the phase diagram corresponding to liquid, 
glassy and nearly crystalline states are mapped out, and the 
nature of the transitions is determined. The liquid to glass 
transition changes from first to second order as the strength 
of the disorder is increased. For weak disorder, the system un- 
dergoes a first order crystallization transition as the density 
is increased. Beyond a critical value of the disorder strength, 
this transition is replaced by a continuous glass transition. 
Our numerical results are compared with those of analytical 
work on the same system. Implications of our results for the 
field-temperature phase diagram of type-II superconductors 
are discussed. 
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I. INTRODUCTION 

The equilibrium phase diagram of a classical system 
of interacting particles in a quenched, random, pinning 
potential is an important subject on which much effort is 
currently being spent [Q. There are several experimen- 
tally studied systems, such as vortices in the mixed phase 
of high-T c superconductors , fluids confined in porous 
media ||, magnetic bubble arrays M, and Wigner crys- 
tals H , which provide physical realizations of a collection 
of interacting classical objects in the presence of an exter- 
nal, time-independent, random potential. In the absence 
of such a potential, systems of this kind are expected 
to crystallize at sufficiently low temperatures. Several 
years ago, Larkin g showed that the presence of arbi- 
trarily small amount of random pinning disorder destroys 
long-range translational order in all physical dimensions 
d < 4. However, recent theoretical studies sug- 
gest that weak disorder distorts the crystalline state only 
slightly, leading to a phase with perfect topological order 
and logarithmic fluctuations of the relevant displacement 
field. This phase, with quasi-long-range translational or- 
der and power-law Bragg peaks in the structure factor, is 



called a "Bragg glass" ||. The transition point between 
a Bragg glass and the high-temperature liquid phase is 
likely to be shifted with increasing disorder, but the tran- 
sition is believed to remain first order as long as the dis- 
order is weak. A question of obvious interest is how this 
transition temperature and the nature of the transition 
depend on the strength of the random potential. 

As the relative strength of the disorder is increased, 
so that the week-disorder situation described above no 
longer applies, the Bragg glass phase is expected to un- 
dergo a transition to a topologically disordered amor- 
phous phase with only short-range translational correla- 
tions. It is not yet clear whether this phase is thermody- 
namically distinct from the high-temperature liquid. An 
interesting possibility is that it is analogous to the glassy 
phase obtained by supercooling a liquid below the struc- 
tural glass transition temperature in the absence of ex- 
ternal quenched disorder ||. If this is so, then the phase 
diagram of such systems would contain three phases: a 
Bragg glass phase obtained at low temperature and weak 
disorder, an amorphous (without quasi-long-range trans- 
lational order) glassy phase at low temperatures and 
strong disorder, and a weakly inhomogeneous (because 
of the presence of the external random potential) liquid 
phase at high temperatures. The glassy phase would be 
a thermodynamically stable one in these systems. This 
is different from the situation in the absence of external 
disorder where the crystalline state is the true equilib- 
rium state near the structural glass transition and both 
the supercooled liquid and the glass are metastable. In 
other words, the presence of external disorder may lead 
to the possibility of occurrence of a true, thermodynam- 
ically stable, glassy phase. 

The phase diagram in the temperature (T) - mag- 
netic field (H) plane of layered, highly anisotropic, type- 
II superconductors such as Bi2Sr2CaCu20s in a magnetic 
field perpendicular to the layers is a credible candidate 
to exhibit these three phases. For a wide range of values 
of H, the flux lines in these materials may be regarded as 
columns of interacting "pancake" vortices residing on 
the layers, and the properties of the mixed phase may be 
described in terms of the classical statistical mechanics 
of these point-like objects. In these compounds, at low 
enough fields, a flux-lattice melting transition separates 
a nearly crystalline state of the flux lines from a disor- 
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dered "vortex liquid" state. The first-order character of 
this transition has been carefully documented [jl0| . When 
the magnetic field is increased, the transition becomes 
continuous jnj[n], and the nearly crystalline state ap- 
pears to be replaced by an amorphous state called "vor- 
tex glass" which is endowed with glassy properties 
such as non-ohmic current-voltage characteristics jl3| . It 
is generally assumed |l2| that the vortex glass phase owes 
its existence to the presence of point-like pinning disor- 
der. Observation of Bragg peaks in neutron scattering 
experiments [Q confirms that the phase at low H and T 
is a Bragg glass. As the effective strength of the disorder 
is increased, either indirectly by increasing H (which is 
believed to increase [|ll| the effective strength of the dis- 
order), or directly by increasing the amount of defects in 
the sample , the Bragg glass phase changes over to the 
vortex glass. The latter is separated from the liquid by 
a continuous transition ]l6[ |. This phase diagram, thus, 
suggests that the first-order liquid-to-crystal transition 
in a three-dimensional (3d) system of point-like objects 
may be driven by the pinning disorder into a continuous 
liquid-to-glass transition. 

The formation of a glassy phase at strong disorder 
was investigated recently |l7j in an analytic study of the 
phase diagram of a system of hard spheres in a random 
pinning potential of arbitrary strength. This work used a 
combination of two "mean-field"- type approaches based 
on the "replicated liquid formalism" ||,|l8],[l9) : the replica 
method [^0| was used for treating the effects of quenched 
disorder, and the hypernetted chain approximation |2l| l 
to calculate the equilibrium correlation functions in the 
liquid in the presence of the pinning potential. These cor- 
relation functions were then the input in a replicated den- 
sity functional jl8| of the Ramakrishnan-Yussouff (RY) 
form [^2| from which the location of the freezing transi- 
tion of the liquid into a nearly crystalline (Bragg glass) 
phase was obtained. The possibility of a liquid-to-glass 
transition was investigated using the phenomenological 
approach of Mezard and Parisi |19|| . The resulting JIt]] 
phase diagram in the density - disorder plane (the den- 
sity, rather than the temperature, is the appropriate con- 
trol parameter for a hard-sphere system) shows three 
thermodynamic phases: a nearly crystalline Bragg glass, 
an amorphous glassy phase, and a low-density liquid. It 
is consistent with the expectation (from earlier work |Ig| ] 
and the Lindemann criterion ]23[]) that the density at 
which the Bragg glass to liquid transition occurs should 
move to higher values as the strength of the disorder 
is increased. The first-order crystallization transition is 
replaced by a continuous glass transition as the disorder 
strength is increased above a threshold value. This phase 
diagram is thus qualitatively similar to that proposed for 
some layered type-II superconductors if, as noted above, 
the density is replaced by the temperature T and the 
disorder strength by the magnetic field H. 

In the present paper, we report the results of a numeri- 
cal investigation of the phase diagram of the same system, 
i.e. a hard-sphere fluid in the presence of a random pin- 



ning potential with short-range spatial correlations. Our 
work involves the use of direct numerical minimization to 
study the effects of the presence of an external random 
potential on the minima of a discretized version of the RY 
free energy functional for the hard-sphere system. It is 
known |2J] that in the absence of external disorder, this 
model free energy functional exhibits, at sufficienly high 
densities, a large number of "glassy" local minima char- 
acterized by inhomogeneous but aperiodic density distri- 
butions. In addition, a global minimum corresponding 
to the crystalline solid is also found at high densities if 
the sample size and the discretization scale are commen- 
surate with the crystal structure. If they are incommen- 
surate, only the glassy minima are present. We have 
carried out extensive numerical investigations of the re- 
sulting free-energy landscape p5|-p8| in the absence of 
disorder. In the present study, and with the physical 
situations described above in mind, we develop similar 
numerical methods to find the location and structure of 
the local minima of the same model free energy with the 
addition of the presence of a time- independent, random, 
one-body potential. 

Using these numerical methods we investigate how the 
uniform liquid, crystalline solid and glassy minima of the 
free energy in the absence of the random potential evolve 
as the strength of the random potential is increased from 
zero. We also examine the dependence of the free en- 
ergy of these minima, and of the density structure (as 
given by the two-point density correlation function) of 
the system at these minima, on the strength of the dis- 
order. In this picture, a transition from one phase to 
another is signalled by the crossing of the free energies of 
the corresponding minima of the free energy. By moni- 
toring where these crossings occur as the density and the 
strength of the disorder are varied, we are able to map 
out the phase diagram in the density - disorder plane. 
This phase diagram is qualitatively very similar to the 
one obtained in the analytic work. For weak disorder 
we find, in the commensurate case as described above 
where a crystalline minimum exists, a first-order liquid- 
to-crystal transition that moves to higher density as the 
strength of the disorder is increased. In the metastable 
"supercompressed" regime (i.e. at a density higher than 
the value at which equilibrium crystallization takes place 
for the commensurate case), we find in all cases a liquid- 
to-glass transition. The density at which this transition 
occurs decreases (very slowly for the largest systems stud- 
ied, which are incommensurate, and more rapidly for the 
smaller, commensurate systems) as the strength of the 
disorder is increased. The nature of this glass transi- 
tion depends on the strength of the disorder: it is first 
order when the disorder is weak, but it changes to sec- 
ond order beyond a certain critical value of the disorder 
strength. For the commensurate case, the crystallization 
line crosses the glass transition line at or very near the 
same critical value of the disorder strength, so that the 
system at stronger disorder then undergoes a liquid-to- 
glass transition (instead of the liquid-to-crystal transition 
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found for weak disorder) as the density is increased from 
a low value. The continuous nature of the glass transition 
in the large disorder regime is in contrast with the first 
order transition from the liquid to a crystalline or glassy 
state (depending on the commensurability) at small val- 
ues of the disorder strength. Thus, this work provides 
support to the prediction that the first-order liquid-to- 
crystal (Bragg glass) transition should change over to a 
continuous liquid-to-glass transition as the strength of 
the pinning disorder is increased above a critical value. 

The rest of the paper is organized as follows. In sec- 
tion II, we define the model studied here and outline the 
numerical procedure used in this study. The numerical 
results obtained for the different transition lines in the 
density-disorder plane are described in detail in section 
III. Section IV contains a summary of our main results 
and a few concluding remarks. 



II. METHODS 

A. The Free Energy functional 

As discussed in the Introduction, our starting point is 
the free energy as a functional of the time-averaged local 
density p(r) at each point r. We write this free energy 
functional in the form: 



F[p]=F RY [p}+F s [p} 



(1) 



where the first term in the right-hand side is the RY free 
enrgy functional |2^] for hard spheres in the absence of 
disorder, and the second is the contribution arising from 
the presence of a quenched random potential. Thus we 
have: 



PF RY [p] = J dv{p{v) ln(p(r)/ Po ) - 6p(r)} 

-\jdrj dv'C{\v-r'\)5p{v)5p{v'). (2) 

Here, we have defined Sp(r) = p(r) — po as the deviation 
of p(r) from po, the density of the uniform liquid, and 
taken the zero of the free energy at its uniform liquid 
value. In Eq.(||), f3 — l/(/csT), T is the temperature 
and the function C (r) is the direct pair correlation func- 
tion [ pl| of the uniform liquid at density po, which can be 
analytically expressed in terms of the usual dimension- 
less density for hard spheres of diameter a, n* = po<J 3 , 
by making use of the Percus-Yevick approximation 
This approximation is sufficiently accurate in the density 
ranges (n* < 1.0) considered in this paper. We write 
also: 



I3F S [ P ] = / dv8p(v)V s (v) 



(3) 



where V s {r) is an external potential (in dimensionless 
form) representing the random, quenched disorder. We 



will assume that V s has zero mean and short-range Gaus- 
sian correlations as detailed below. 

In order to carry out numerical work, we discretize our 
system. We introduce for this purpose a simple cubic 
computational mesh of size L 3 with periodic boundary 
conditions. On the sites of this mesh, we define density 
variables pi = p(ri)h 3 , where p(vi) is the density at site 
i and h the spacing of the computational mesh. It is 
known from previous work 
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that in the absence of 
any random potential, this discretized system crystallizes 
at sufficiently high densities if the quantities h and L are 
commensurate with a fee structure with appropriate lat- 
tice spacing, whereas no crystalline state exists when the 
computational mesh is incommensurate with a fee struc- 
ture. Both commensurate and incommensurate systems 
exhibit ]24|-p7| many glassy (inhomogeneous but aperi- 
odic) minima of the free energy at densities higher than 
the value at which crystallization occurs in commensu- 
rate samples. 

To model the random potential V s (r), we introduce 
random variables {Vi} defined at the sites of the compu- 
tational mesh. These variables are uncorrelated with one 
another, and distributed according to a Gaussian proba- 
bility distribution with zero mean and variance s. Thus, 
s represents the dimensionless strength of the disorder. 
In terms of these quantities, the dimensionless free energy 
of our discretized system has the form 

PF = ^2{Pi \n(pi/p e ) - (pi - pe)} 

i 

- \ C ij(Pi ~ Pt)(Pj ~ Pi) +X] Vi ( pi ~ Pe )> ^ 

i j i 

where the sums are over all the sites of the computa- 
tional mesh, pi = p$h 3 1 and CV, is the discretized form of 
the direct pair correlation function C (r) of the uniform 
liquid. 

Our objective is to study the phase diagram of this sys- 
tem in the (n* , s) plane, in which in principle crystalline, 
liquid, and glassy phases may be found. The thermody- 
namics of hard spheres in the clean limit is determined 
by the dimensionless density n* only. Our rcscaling of 
the potential V s by f3 (see Eq.(||)) ensures that s is now 
the only additional relevant variable. We wish to locate 
various transition lines in this (n*,s) plane, that is, we 
wish to determine which phase (crystal, glass or liquid) 
is the thermodynamically stable one at different points 
in this plane. We also wish to know when and how a cer- 
tain phase becomes metastable or unstable as we move 
around in the (n* , s) plane. In our mean- field descrip- 
tion, different phases are represented by different minima 
of the free energy. If several local minima of the free en- 
ergy are simultaneously present, then the minimum with 
the lowest free energy represents the thermodynamically 
stable phase and the other local minima correspond to 
metastable phases. A crossing of the free energies of 
two different minima signals a first-order phase transi- 
tion. The point where a minimum becomes locally un- 
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stable (i.e. changes from a true minimum to a saddle 
point or disappears altogether) corresponds to a mean- 
field spinodal point representing the limit of metastabil- 
ity of the corresponding phase. A merging of the tran- 
sition point with the spinodal points of the two phases 
signals a continuous phase transition in this description. 
Thus, a study of how the minima of the free energy of 
Eq.(^) evolve as n* and s are changed is sufficient for 
mapping out the mean-field phase diagram in the (n* , s) 
plane. 

We locate the minima of the free energy by using a 
numerical procedure generalized from that originally de- 
veloped for the clean case [Q . This procedure works by 
changing the local density variables {pi} in a way that 
ensures that these changes always decrease the free en- 
ergy. Given an initial configuration of the variables 
this procedure finds, by constantly moving downhill on 
the free-energy surface in the multidimensional configu- 
ration space spanned by the L 3 variables {pi\, the local 
minimum whose basin of attraction contains the intial 
state. Thus, different local minima of the free energy 
can be located by using this minimization procedure for 
different, appropriately chosen, initial configurations. 

As noted earlier, there are in our system three differ- 
ent kinds of free-energy minima: liquid, crystalline and 
glassy. In the clean limit (s = 0), it is easy to distinguish 
among them: the liquid minimum has uniform density 
(f i = Pi for all i), the crystalline minimum has a periodic 
distribution of the density variables (pi is close to unity 
at mesh points corresponding to the sites of a fee lattice, 
and close to zero at all other mesh points), and a glassy 
minimum exhibits a strongly inhomogeneous nonperiodic 
density distribution (some of the p^s are close to unity 
and the others are close to zero). This symmetry-based 
distinction among minima of different kind becomes less 
clear when the external random potential is turned on: 
for s =/= 0, the density distribution in the liquid phase is 
not completely homogeneous, and the crystalline state is 
not strictly periodic. 

We use here, therefore, a procedure of "adiabatic con- 
tinuation" to distinguish among the liquid, crystalline 
and glassy minima in the presence of the disorder. This 
procedure works as follows: We start with a minimum of 
a particular kind obtained at s = for a given value of 
n* . There is no difficulty in generating the liquid (and 
if appropriate the crystalline) configuration for the pure 
system. Glassy states at s = are easily obtainable also, 
in the right density ranges, by the procedures described 
in Ref. | f27f . Indeed, we have used in many cases the same 
density configurations obtained there, which were avail- 
able as computer files. After thus choosing the initial 
state, we generate a set of uncorrelated random num- 
bers r i7 i = 1, . . . , L 3 , distributed according to a Gaus- 
sian with unit variance. A "realization" of the random 
potential {Vi} is obtained by multiplying these random 
numbers by the strength parameter s. The initial s = 
minimum is then "followed" to finite s by increasing s in 
small steps 5s [E9J . After each step increase, the min- 



imization routine is run, to find the nearest local mini- 
mum. Thus, in the first step of this process, the initial 
configuration is that of the minimum obtained at s = 
and the values of V, are set at (Ss)ri. The resulting den- 
sity configuration obtained from the minimization rou- 
tine is then used as the initial configuration for the next 
step, with the values of V, incremented to 2(5s)ri. Dur- 
ing this process, the random variables {r{\ are held fixed 
- only the strength parameter s in increased in steps of 
5s. By iterating this procedure, minima of different kinds 
obtained at s — for a certain n* are "followed" at con- 
stant density to the desired value of s. We use the terms 
"liquid" , "crystalline" and "glassy" to denote the contin- 
ued s minima obtained from as = minimum of the 
corresponding kind by using this continuation procedure 
without crossing transition lines. We will see below that 
even at large s, there are distinguishable differences in 
the structure of the different kinds of minima. 

Once a minimum of the desired kind is obtained at 
a particular point in the (n*,s) plane, the free energy 
at the minimum reached, as well as the entire density 
configuration of the system at the minimum are obtained 
and can be analyzed. The translational correlations can 
be quantified by the two-point correlation function g(r) 
of the density variables {pi}. This function is defined as 



9(r) 



(5) 



<>./ 



where the distance r is measured in units of a, p = 
J2i Pi/L 3 is the average value of the pi variables at the 
minimum under consideration, and fij(r) = I if the sep- 
aration between mesh points i and j lies between r and 
r + Ar (Ar is a suitably chosen bin size), and fij(r) = 
otherwise. This function represents the spatial correla- 
tion of the time- averaged local density, and is distinct 
from the equal-time, two-point density correlation func- 
tion which is often called g{r) in the literature. We also 
calculate p m ax , the maximum value of the pi variables at 
the minimum, which gives additional information about 
the inhomogeneity when contrasted with p or its rescaled 
equivalent p av = p(a/h) 3 at the minimum. 

In addition to examining the transitions by looking at 
discontinuities in F, g(r) and the density configurations, 
we also directly check on the stability of the correspond- 
ing minima. The stability of a local minimum requires 
that all the eigenvalues of the Hessian matrix M whose 
elements are given by 



Mi 



_ d 2 {/3F) 1 



dpidpj 



— 5i 



Ci 



(6) 



evaluated at the minimum must be positive. This matrix 
is difficult to handle numerically if the minimum under 
consideration is strongly inhomogeneous, with some of 
the pi's very close to zero. In such cases, the l/pi in 
the first term on the right-hand side of Eq.@ causes 
numerical difficulties. To avoid this problem, we consider 
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instead the closely related matrix M' whose elements are 
given by 

M'ij = \ JkM,, n p] = Sij - Cijy/pipj, (7) 

evaluated at the minimum under consideration. It is easy 
to show that an instability of the minimum corresponds 
to the vanishing of the smallest eigenvalue A of this ma- 
trix. In our numerical work, we calculate the value of 
A in order to check whether the minimum under study 
becomes unstable as n* or s is varied. 

In our computations we have included the density 
range from n* = 0.65 to n* = 0.95, and values of s 
from zero to about two. These are sufficient to encom- 
pass the phenomena that we wish to study. We have used 
three lattice sizes, L = 12,15, and 25. For the last two we 
have used an incommensurate ratio h/a — 1/4.6, whereas 
for the smallest lattice we have taken the commensurate 
value h/a = 0.25. 

III. RESULTS 
A. General considerations: Phase diagram 



transition. In the commensurate case, the density n* c 
is above the crystallization density n* D , and the free en- 
ergy of the crystalline minimum is lower than that of the 
glassy minima. Thus, the glass transition in the pure 
system occurs in a "supercompressed" regime where the 
crystalline state is the thermodynamically stable one. 

When we include the effects of the disorder (s > 0), 
we find yet another density, n* A , at which the liquid min- 
imum becomes locally unstable (i.e. ceases to exist as a 
local minimum of the free energy). For weak disorder, 
the value of n* A is large (substantially higher than n* B ) 
so that the four densities n* A , n* B , n* c , and n* D are in 
decreasing order. Thus, we have four (three in the in- 
commensurate case where the crystalline state is absent) 
functions n^(s) with X = A, B,C, D, representing pre- 
cisely the four transitions or instabilities defined above. 
We denote the corresponding lines in the (n*, s) plane as 
the A, B, C, D lines. The determination of the location of 
these lines is one the main results of our work. These re- 
sults will be discussed below, but to fix ideas and to make 
this discussion easier to follow, we show in Fig. [j] these 
four lines for the L = 12 commensurate case. There, 



At a general point in the (n*, s) plane, the system may 
have a number of local minima, one of which is the ab- 
solute minimum while the others are metastable. Possi- 
ble minima are that corresponding to the liquid (the one 
with uniform density at s = and its continuation to 
finite s), the crystalline minimum, by which we similarly 
mean the one with a periodic structure at s = 0, and its 
continuation as described in the preceding section, and 
a large number of glassy minima. As the values of n* 
or s change, free energy minima may in general appear 
and disappear, and the free energy values of those that 
remain change. There will therefore be a number of in- 
stabilities and transitions, which are the main subject of 
our study. 

Consider first the previously studied |24 2?]j2^] case of 
the disorder- free system (s = line). There, only the 
uniform liquid minimum is present at low densities. As 
n* increases, a crystalline minimum appears if the com- 
putational mesh is commensurate. When n* is further 
increased, a density is reached at which the crystal be- 
comes thermodynamically stable, that is, its free energy 
becomes lower than that of the liquid state. We will de- 
note this density as n * D . Regardless of commensurability, 
many glassy minima appear as the density is further in- 
creased. We denote by n* c the density at which the first 
glassy minimum makes its appearance. Alternatively, 
one may consider the evolution of the glassy minima as 
n* is decreased from a large initial value, and define n* c 
as the density at which the last remaining glassy mini- 
mum becomes locally unstable and disappears: the free 
energy of this last remaining glassy minimum crosses that 
of the liquid at a density n* B which is somewhat higher 




than 



This density corresponds to a liquid to glass 



FIG. 1. The overall phase diagram of the hard sphere sys- 
tem in the density (n*) - disorder (s) plane, obtained for the 
L — 12 commensurate sample. The meaning of the line labels 
is explained in the text. The results shown are averages over 
5 realizations of the disorder. 

the general structure of the phase diagram, includ- 
ing the general shape of the four lines n^(s) can be 
seen. Similarly, we show in Fig. ^| the three lines 
n* x {s),X ~ A,B,C (from top to bottom) found in the 
incommesurate, L = 25 system. The similarities and dif- 
ferences between the comensurate and incommensurate 
cases are discussed below. The lines in the phase diagram 
for the incommensurate L — 15 case are within error bars 
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the same as those shown in Fig. |[ so that the differences 
between Figs. [j] and || must be attributed to different 
commensurability rather than to different sample size. 
There are certain trends that can be easily discerned 
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FIG. 2. The overall phase diagram for the incommensurate 
case at size L = 25 as explained in the text. The diamonds 
represent the A line, the crosses the B line and the dashed 
line is the C line. Sample error bars have been indicated. 
They reflect sample-to-sample variations for six to twelve (the 
number increases with s) realizations of the disorder. 
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FIG. 3. Liquid phase correlations. The pair correlation 
function g{r) as defined in Eq. (|H]) plotted as a function of 
r, the distance in units of a, for the liquid- like minimum at 
density n* = 0.66. The curves shown, in order of increasing 
peak height at r = 1, correspond to s = 0.2,0.6, 1.0, 1.4, 1.8. 
The system size is L — 25. 



when one follows a free energy minimum as s is increased 
at constant n*. If one starts from the uniform liquid 
minimum at s = and a relatively small value of n* , the 
free energy value at the minimum (initially zero according 
to our convention) decreases steadily with increasing s. 
For the case in Fig. || at n* = 0.66, for example, (3F 
is close to —180 at s = 1.8. The density distribution 
becomes progressively less uniform, with p m ax, which at 
s = equals the average value p = pt, rising by more 
than one order of magnitude as s increases from zero to 
one. For a deep glassy state at a relatively large value of 
n* , the free energy is strongly negative even at s — 0, and 
its value decreases further as s increases. For example, at 
n* = 0.78 the glassy minimum for L = 25 with (3F = —63 
at s = can be|QO^gtinued to a minimum with (3F equal 
to —231 at s = 1.8. The density distribution at a glassy 
minimum is considerably more inhomogeneous than that 
of the liquid minimum continued to the same value of s 
and it is less sensitive to the value of s: the quenched 
disorder has less effect on a state that is inhomogeneous 
and disordered to begin with. 

These trends in the behavior of liquid and glassy min- 
ima as s is increased from zero are clearly illustrated by 
examining the pair correlation function g(r), defined in 




0.5 1 1.5 2 2.5 3 3.5 4 4.5 
r 

FIG. 4. The pair correlation function g(r) for a glassy min- 
imum. The curves shown correspond to the same values of L 
and s as in Fig.^, but for a glassy minimum at n* = 0.78 as 
discussed in the text. 

Eq.([5]), at each minimum. In Fig. ^ we show g(r) com- 
puted for the liquid minimum at size L = 25 and density 
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n* = 0.66. The curves shown, in order of increasing value 
of the peak near r = 1, correspond to increasing values 
of s = 0.2,0.6, 1.0, 1.4, 1.8. There is a clearly visible in- 
crease in structure, which becomes more evident as the 
value of s increases beyond unity. However, this level of 
structure is still quantitatively different from that found 
for glassy minima at relatively high densities. This can 
be seen by comparing Fig. |^ with Fig. |J where we plot 
g(r) for a L = 25 glassy minimum continued from s = 
to s = 1.8 at n* = 0.78, as mentioned in the preceding 
paragraph. We see that the s-dependence of the struc- 
ture is now much weaker, and the heights of the peaks 
at r = and near r = 1 are much larger than those in 
Fig. ||. These results can be compared to those found in 
the replica calculation To make contact with those 
results, our g(r) for the liquid-like minimum should be 
compared with the function g$ (r) of the replica symmet- 
ric solution, and our g(r) for a glassy minimum with the 
function gi(r) of the replica-symmetry-broken solution. 
Although, due to differences in the modeling of the ran- 
dom potential and effects of discretization in the present 



study (some of these effects are discussed in section IV), 
a detailed, quantitative comparison of our results with 
those of Ref. |l7j is not possible, it is clear that the main 
features we have discussed are qualitatively similar. 

The crystalline minimum obtained for s = in com- 
mensurate systems at sufficiently high densities shows 
very little change in structure as it is followed to non- 
zero values of s. Any effects of weak pinning disorder on 
the crystalline order may be too subtle [Ug] to show up 
at the system sizes and discretization scales used here in 
the commensurate case. 



B. Instability of the liquid minimum 

We consider first the A line, that is, the density at 
which the liquid minimum becomes locally unstable as 
n* is increased from a low initial value, keeping s fixed. 
This transition is detected at any desired value of s in the 
following way. At a density previously determined to be 
well below the value of n* A (s) (this determination is eas- 
ily performed by trial and error) , one "follows" the s = 
liquid minimum, as previously explained, to the value of 
the disorder strength being studied. The density configu- 
ration at this minimum is the initial condition. Then, one 
proceeds to increase n* by small intervals, thus moving 
up along a vertical line in the phase diagram. At every 
value of n* that is reached, we run our minimization rou- 
tine (using the configuration at the minimum obtained 
at the previous step as the starting point) to locate the 
nearest minimum. The density configuration at the min- 
imum is analyzed and then used as the initial condition 
to study the next higher density. 

In the initial stages of this process, the system remains 
in the liquid-like minimum, with little change in its prop- 
erties. However, as n* reaches the value n* A (s), discon- 



tinuities are found. These are more prominent for the 
larger system sizes (Fig. ||) and particularly dramatic for 
values of s not too large. The free energy drops abruptly, 
as the liquid minimum becomes unstable and the system 
has to find some other nearby minimum (our numerical 
minimization procedure is designed to converge only to 
stable local minima of the free energy) . Computationally, 
this is heralded by a very sharp and obvious increase in 
the number of iterations required by our numerical pro- 
cedure to find the free energy minimum nearest to the 
starting configuration. This new minimum is invariably 
glassy, as one might expect, since a considerable number 
of glassy minima are close in configuration space to the 
liquid- like minimum |^7j . The value of the free energy at 
the minimum that the system has reached drops sharply 
as the n A (s) value is crossed, because the free energies 
of glassy minima are considerably lower in the region of 
the (n* , s) plane being considered. Also, every measure 
of structure in the system inceases abruptly, since, as 
discussed above in connection with Figs. || and f|, glassy 
states are much more inhomogeneous than the liquid-like 
ones in this region of the (n* , s) plane. 

An example of the behavior found is displayed in 
Figs. |5| and [|. In the main part of Fig. pi we show 
the evolution of the free energy as n* is increased 
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FIG. 5. Discontinuities at the A line. In the main plot, 
the free energy in dimensionless form for a L = 12 sample 
(s = 0.6) is plotted as a function of density. A sharp drop in 
the free energy is seen as the liquid minimum becomes unsta- 
ble and the system switches to a glassy minimum. As shown 
in the inset, this switch is also reflected in the discontinuity in 
A, the smallest eigenvalue of the matrix M' defined in Eq.(Q). 

in steps of 0.001, keeping s fixed at 0.6 for a L = 12 
sample. One can clearly see that f3F varies little while 
the system remains in the liquid minimum and jumps 
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abruptly as this minimum becomes unstable near n* A ~ 
0.78. The behavior for the larger incommensurate sam- 
ples is quite similar, the main difference being that the 
drop in (3F is much larger, and that the transition oc- 
curs, at this value of s, at n A ~ 0.84 for both L = 15 
and L = 25. The value of n* A can readily be found to 
very high precision and it varies little as one averages 
over different realizations of the quenched disorder, for 
the same s. The error bars shown in Fig. |^ correspond to 
an average over six to twelve realizations (the larger num- 
ber at larger s.) The results in Fig. [I] are averages over 
five realizations. In the inset, we show that the smallest 
eigenvalue, A, of the matrix M' defined in Eq.(^) ap- 
proaches zero as n* approaches n* A from below. This is 
as would be expected - as noted in section [o], the instabil- 
ity of a local minimum is signalled by the vanishing of A. 



short-range order present in a glassy minimum. This can 
also be seen from Figs. || and f|. The value of p max also 
increases by a considerable amount, indicating the in- 
creased inhomogeneity of a glassy minimum relative to 
the liquid-like one. The small increase in the value of p av 
reflects that the average density at a glassy minimum is 
slightly higher than that at the liquid-like minimum. 

The behavior discussed above changes as s is increased. 
The change occurs near s — 1 for L — 12, and at some- 
what larger s for the other system sizes, as the A, B, C 
lines become very close. The results obtained in the 
largcr-s regime are described in the next subsection. 

C. Instability of glassy minima and the liquid to 
glass transition 
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FIG. 6. Example of how the system becomes more struc- 
tured as the A line is crossed for an L = 12 sample at s = 0.6. 
The height g m ax of the first finite-r peak in g(r) increases dis- 
continuously, the density nonuniformity represented by p ma x 
exhibits a large increase, and the average density p av shows 
a small discontinuous increase. In order to be able to use a 
single vertical scale, we have displayed p av rather than p (see 
text). 

In Fig. |(| we show three quantities which characterize the 
nature of the density distribution at a minimum. These 
are: g m ax, the value of the pair correlation function g(r) 
at its first finite r maximum (near r = 1); p m ax, the 
maximum value of the pi] and the dimcnsionlcss aver- 
age density p av defined in section All these quanti- 
ties exhibit discontinuous changes as the system switches 
minima at n* — n* A ~ 0.78, (or at n* A ~ 0.84 in the in- 
commensurate case). g m ax remains close to unity as long 
as the system stays in the liquid state, and then jumps to 
a substantially larger value consistent with the increased 



To find the B and C lines, we must start with a care- 
fully chosen glassy configuration at a relatively high n* 
and fixed s, and then follow this configuration to lower 
densities by decreasing n* in small steps (5n* ~ 0.001), 
keeping the value of s unchanged. This is continued un- 
til the minimum becomes unstable and the minimization 
routine converges to a new minimum which, if the start- 
ing minimum is chosen as described below, turns out to 
be the liquid-like one. The density at which this occurs 
defines the value of n* c . Then, comparing the free energy 
of the glassy minimum with that of the liquid minimum 
obtained for the same realization of the disorder, it is 
easy to determine the value of n* B - this is the value of 
n* at which the two free energies are equal. 

The determination of the appropriate starting glassy 
minimum is nontrivial. Glassy minima for s ^ are 
obtained by continuation from those of the pure system 
(s = 0). One may think that the best choice would be 
to take the glassy minimum with the lowest free energy 
at the starting (n*, s) point. In practice, this is difficult 
to implement because an exhaustive enumeration of all 
the glassy minima is computationally very hard. The 
glassy minimum with the lowest free energy at a particu- 
lar point in the (n* , s) plane does not in general continue 
to have the lowest free energy as the values of n* and s are 
changed. Also, in the pure system, all the configurations 
obtained by applying one of the symmetry operations of 
the computational mesh to the density configuration at a 
particular glassy minimum also correspond to local min- 
ima with exactly the same free energy. For s ^ 0, all 
these symmetry-related minima have to be considered 
separately because the presence of the random potential 
destroys the symmetries present in the pure limit. 

We have not found a rigorous solution to this prob- 
lem. Instead, we first carried out an exploratory study of 
how the locations of the B and C lines in the phase dia- 
gram depend on the choice of the initial glassy minimum. 
The following choices, among others, were considered in 
our initial exploration, (a) One of the low-lying s = 
glassy minima, continued to finite s. (b) Beginning with 
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the same starting configuration as in (a) and a specific 
realization of the random variables {Vi}, minimize the 
random potential energy (the last term in Eq.(Q)) with 
respect to all symmetry operations of the computational 
mesh. This attempts to find the configuration that min- 
imizes, among all the symmetry related ones, the con- 
tribution of the random potential to the free energy but 
it is not quite rigorous because the minimization is per- 
formed using the values of {pi} at the s — minimum. 
(c) The glassy minima to which the system moves when 
the density is increased above the A line, as discussed in 
the preceding subsection. 

The outcome of this study is that the locations of the 
B and C lines in the (n* , s) plane are not sensitive to the 
choice of the glassy minimum as long as it is one of the 
low-lying minima. (Even when we have deliberately or 
accidentally chosen a "wrong" , non-low-lying, minimum, 
we have found that the system often spontaneously makes 
a glass-to-glass transition p0| to a low-lying minimum as 
one decreases n* above the B line.) The variation of the 
values of n* B and n* c for different choices of the glassy 
minimum is comparable to the uncertainty of these val- 
ues arising from sample-to-sample variations caused by 
differences in the realization of the disorder. The re- 
sults described below were obtained (unless otherwise in- 
dicated) from runs in which a low-lying glassy minimum 
obtained from continuation of one at s = was taken 
to be the initial state for the density-lowering run. As 
explained at the beginning of this subsection, the initial 
configuration is followed to lower densities at fixed s and 
the n^(s) and ng(s) points are found for that value of s. 
For relatively small values of s, the signatures of the C 
instability are very easy to detect: they are similar to the 
discontinuities shown in Figs. |5| and [j[ At larger values 
of s more care is required. 

For small values of s, the A, B and C lines are well 
separated from one another. However, as the value of 
s is increased, these three lines begin to approach each 
other. As shown in Figs. [I] and ^, the separation between 
lines A and B decreases rather rapidly with increasing 
s, while the separation between lines B and C decreases 
more slowly. Finally, near s = 1, these three lines appear 
to merge with one another for the L = 12 system. For 
L = 25 (and also for L = 15), the separation between 
them does not exceed the combined error bars, but sep- 
arate B and C transitions can be detected in most (not 
all) "runs" (i.e. realizations of the disorder) as explained 
in detail below. At larger s, it becomes increasingly dif- 
ficult to resolve these three lines as they come close to 
each other. Since lines A and C represent, respectively, 
the limits of stability of the liquid and glassy minima and 
line B corresponds to the first-order liquid-glass transi- 
tion, a merging of these three lines suggests that this 
transition becomes continuous as s is increased beyond a 
"tricritical" value which would be close to unity for the 
L = 12 commensurate sample and somewhat larger for 
the incommensurate samples. Another possibility is that 
the first-order liquid-glass transition disappears beyond 



a critical point near s = 1. 

To examine the behavior in this region more closely, we 
have carried out several numerical experiments in which 
the value of n* is "cycled" through the liquid-glass transi- 
tion, keeping s fixed at values close to unity. In this way, 
the three lines are detected in the same "run" . These 
numerical experiments are similar to simulations of hys- 
teresis in magnetic phase transitions. We start with the 
liquid minimum at a low value of n* (below line C), and 
increase n* in small steps, keeping s fixed. The liquid 
minimum is thus followed to higher densities until it un- 
dergoes a rapid change signalling a possible instability. 
The process of increasing n* in small increments is con- 
tinued for a few more steps, and then the local minimum 
so obtained is followed to lower densities by decreasing 
n* is small steps. This is continued until the starting 
value of n* is reached. If the liquid-glass transition at 
the chosen value of s is first-order with the three densi- 
ties n* Al n* B , and n* c separated from one another, then the 
cycling experiment described above should exhibit clear 
evidence of hysteresis. This is indeed what we find, for 
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FIG. 7. Hysteresis and discontinuities across the liq- 
uid-glass transition at small values of s. In the main plot, 
the dimensionless free energy of the stable minimum is plot- 
ted vs. density as one cycles across the A, B and C lines as 
explained in the text. Hysteresis is clearly observed. In the 
inset, the quantity g max is shown. The results shown are at 
s = 0.8 for a commensurate L = 12 system, but the same be- 
havior is found in this range of s for incommensurate systems. 

all system sizes and at every run, if the value of s is 
lower than a certain critical value. A typical example is 
shown in Fig. [7] which shows the results for a L = 12 
sample at s = 0.8. The hysteresis in the free energy and 
gmax (shown in the inset) is evident: the liquid mini- 
mum becomes unstable at n* A ~ 0.735 as n* is increased 







from a low initial value, while the glassy minimum found 
for n* > n* A can be continued all the way down to 
n* c ~ 0.720 before it becomes unstable. The liquid-glass 
transition occurs at n* B ~ 0.725 where the two branches 
of the free energy cross. The same situation occurs for 
the incommensurate L = 25 system except that the val- 
ues of the transition points are n* A ~ 0.79, n* B ~ 0.73 
and n* c ~ 0.71 for s = 0.8. The results at L = 15 are, 
within error bars, the same as those for L = 25 at this 
value of s. 

The behavior in Fig. [?] is to be contrasted with that 
shown in Fig. || which displays the results of the cycling 
experiment on a L = 12 sample at s = 1. The distribu- 
tion of the random variables {r{\ in this sample is the 
same as that of Fig - only the strength of the disor- 
der is changed. In this figure, there is no evidence of 
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FIG. 8. Cycling across the liquid-glass transitions for a = 
1 in a L = 12 commensuarte system. The same quantities are 
plotted as in Fig. (?|, and now no hysteresis is seen. Incommen- 
surate systems exhibit the same behavior at somewhat larger 
values of s, but not in all runs. 

hysteresis in the free energy. The plot of g max shown 
in the inset exhibits a sharp change near n* — 0.706 for 
both inreasing-7i* and decreasing-?!* runs, and the results 
for the two runs are nearly identical. Given the rounding 
off errors associated with the numerical procedures we 
use, the small differences between the increasing-n* and 
decreasing-n* values of g m ax are likely to be insignifi- 
cant. We, therefore, conclude that at least within the 
resolution of our numerical procedures, there is no hys- 
teresis at s = 1.0 for this L = 12 sample. This implies 
that the first order transition found in this sample for 
s = 0.8 either becomes a continuous one or disappears as 
the value of s is increased to 1.0. The sharp change in 
the value of g ma x near n* — 0.706 suggests that the tran- 



sition persists as a continuous one. To investigate this 
further, we have calculated the derivatives of g ma x, Pmax, 
and ptot = Pi with respect to n* in the region where 
these quantities change rapidly. We have also exam- 
ined the behavior of A as a function of n* in this region. 
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FIG. 9. Derivatives with respect to n* of the quantities 
Ptot, pmax and g ma x, as defined in the text, plotted as func- 
tions of n* across a putatively continuous liquid-glass transi- 
tion in a L — 12 sample with s — 1.0. The three quantities 
have sharp peaks at n* — 0.706. The eigenvalue A, also de- 
fined in the text, shows a pronounced dip at the same point. 

Results for these quantities are shown in Fig. || for the 
same sample as that of Fig. ||. All the derivatives ex- 
hibit sharp peaks at n* = 0.706, and the value of A goes 
through a minimum that is very close to zero at the same 
point. These results strongly suggest the occurrence of 
a continuous phase transition at n* — 0.706. However, 
due to the limited resolution of our numerical calcula- 
tions and the smallness of sample size, we can not rule 
out the possibility that the observed behavior reflects a 
sharp crossover rather than a true phase transition. Sim- 
ilar results are found for larger values of s. The contin- 
uation of the "transition line" beyond the point where 
the lines A, B and C come together is determined by lo- 
cating the value of n* at which the eigenvalue A reaches 
a minimum. The value of s at which the A, B and C 
lines merge and the hysteresis in the cycling experiment 
disappears is found to be weakly dependent on the real- 
ization of the disorder - it varies between 1.0 and 1.2 for 
the five different L = 12 samples studied. 

For the incommensurate samples, the situation is 
somewhat more ambiguous. For L = 25, the same cycling 
procedure shows that the transition is clearly hysteretic 
for all runs with s < 1.1. For larger values of s, an in- 
creasingly larger percentage of the runs is non-hysteretic 
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(i.e. the results for f3F look like those in Fig. |J), while 
the other runs display a behavior similar to that in Fig. Q 
but with much smaller discontinuities. As s is increased 
beyond s = 1.8, it becomes, in most of the "runs", im- 
possible to distinguish the discontinuities, if any, from 
computer noise. Thus, it is possible in this case to plot 
separate A, B, and C lines all the way up to s = 1.8. 
This accounts for the obvious difference in this respect 
between Figs. |l| and |^. The results for L — 15 are quite 
consistent with those for L = 25, but the smaller sys- 
tem size makes all interpretations more difficult. Thus, 
it is more difficult to identify the precise position of any 
well-defined tricritial point (or a critical point) from the 
results for the incommensurate samples. One might al- 
ternatively say that these incommensurate results are in- 
dicative of a crossover. It is not possible to completely 
rule out that the behavior is different for the commen- 
surate and incommensurate samples, or that the poorer 
resolution of the smaller samples masks discontinuous be- 
havior in some of the larger s runs. 



D. Crystallization 

To study how the crystallization density n* D changes 
as s is increased from zero (i.e. the location of the 
line D), we start with the crystalline minimum ob- 
tained for a commensurate sample at a large value 
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figuration to the desired value of s. This configuration 
is then continued to smaller values of n* by decreasing 
n* in small steps. The crystalline minimum turns out 
to be quite robust under changes of the density and the 
strength of the disorder - the minimization routine con- 
verges rather quickly to the new minimum as the value of 
n* or s is changed by a small amount. While decreasing 
the value of n* , we keep track of the free energy of the 
crystalline minimum and find the value of n* at which 
this free energy crosses that of the liquid minimum for 
the same realization of the disorder. For relatively small 
values of s, the crossing point determines the value of n* D 
for the chosen value of s. Typical results for the cross- 
ing of these two free energies are shown in Fig. [h]. Our 
results for line D, averaged over five realizations of the 
disorder, are shown in Fig. ^. The crystallization transi- 
tion is strongly first order for all values of s. In the small 
s regime, the crystalline minimum has the lowest free en- 
ergy for all densities above line D. Therefore, the lines 
A, B and C do not have any equilibrium thermodynamic 
significance in this regime for a commensurate system: 
the liquid-glass transition at line B can be observed only 
if the crystallization transition at line D is avoided, e.g. 
by rapid compression. 

As shown in Fig. pi the crystallization line crosses the 
liquid-glass transition line at a point which is very close 
to that where the lines A, B and C seem to come to- 
gether. Beyond this point, line D is determined by the 
crossing of the free energies of glassy and crystalline min- 
ima. The procedure is quite analogous to that shown in 
Fig. |l0|. This line, therefore, represents a first order tran- 
sition between crystalline and glassy states in this regime. 
The phase diagram of Fig. [I] implies that the system un- 
dergoes a first order liquid-to-crystal transition for small 
values of s as the density is increased from a low initial 
value. However, as the value of s is increased above a crit- 
ical value (which is close to unity for the L = 12 system), 
the transition as n* is increased becomes a continuous 
liquid-to-glass transition (or perhaps a sharp crossover). 
The glassy state then undergoes a first order transition 
to the crystalline state as the density is increased further. 
The observed curvature of line D for large s also implies 
that the system would undergo a first order crystal-to- 
glass transition as the strength of the disorder is increased 
at constant density. 



IV. SUMMARY AND DISCUSSION 



FIG. 10. Free energy crossing at the crystallization transi- 
tion. The solid and dotted lines represent, respectively, the 
free energies of the crystalline and liquid-like minima of a 
L = 12 sample with s = 0.6. Their crossing point is the 
density n* D (s). 

of n* . We then find the symmetry related configuration 
that minimizes the random potential energy for a par- 
ticular realization of the disorder and continue this con- 



We have mapped out the mean-field phase diagram 
of a hard sphere system in the presence of a quenched 
random potential by numerically studying the evolution 
of the minima of a model free energy as a function of 
the density n* and the strength s of the disorder. The 
phase diagram in the (n* , s) plane exhibits liquid, glassy 
and crystalline (for commensurate samples) phases. We 
find that the standard first order crystallization transi- 
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tion which occurs at s = upon increasing n* retains 
its character at small s as a first order transition from 
a weakly inhomogeneous liquid phase to a nearly crys- 
talline state. The density at which this transition occurs 
increases somewhat with the strength of the disorder. 
We also find for all samples a liquid to glass transition in 
the metastable, "supercompressed" regime. This transi- 
tion is first-order for small s, but within the resolution 
of our results it appears to become continuous as s is 
increased above a critical value. This critical value is 
larger for incommensurate samples. The crystallization 
line in the (n* , s) plane crosses the glass transition line 
near the point where the glass transition becomes con- 
tinuous. Thus, the first order crystallization transition 
found for small s as the density is increased from a small 
initial value is replaced, for sufficiently large s, by a con- 
tinuous liquid to glass transition. The phase diagram also 
shows, at larger s, a first order crystal to glass transition 
as the strength of the disorder is increased at constant 
density. 

All the qualitative features of our phase diagram (i.e. 
its topology, the shapes of the transition and instabil- 
ity lines, and the nature of the transitions) are identical 
to those found in the replica-based analytic study ll7p 
of the same syatem Two of our most important 

results, namely the change in the nature of the liquid to 
glass transition as the strength of the disorder is increased 
above a critical value, and the crossing of the crystalliza- 
tion and glass transition lines above this critical value 
of the disorder strength, were also found in the analytic 
calculation. This similarity between the results of two 
studies using extremely different methodologies strongly 
suggests that the qualitative features of our phase dia- 
gram are correct, at least at the mean-field level. The 
considerable quantitative differences that exist between 
the numerical and replica results, i.e. that all the tran- 
sition and instability lines in our phase diagrams lie at 
substantially lower densities than those obtained in the 
analytic study, have the same origin as the discrepancy 
between our s = results and those of molecular dy- 
namics simulations [fs^| of the pure hard sphere system. 
As noted in our earlier work (2^^], these differences re- 
sult from the discretization of the free energy functional. 
The use of a simple cubic mesh of spacing h ~ 0.2a in the 
discretization procedure increases the relative stability of 
inhomogeneous local minima of the free energy and thus 
leads to substantially lower values for the densities at 
which crystallization and the glass transition occur. The 
quantitative differences between our results in Fig. ^ and 
those in Fig. ||, on the other hand, appear to arise chiefly 
from the incommensurability of the latter sample, rather 
than from the slight difference in the values of h, or even 
from that in the values of L: we have found negligible 
sample-size effects in comparing the L = 15 results to 
those at L = 25 at the same value of h. The effects of 
discretization would presumably disappear for h much 
smaller than the width of the approximately Gaussian 
density distributions near the points where the particles 



are localized at an inhomogeneous minimum of the con- 
tinuum free energy functional. Unfortunately, a numer- 
ical calculation with such small values (~ O.OIct) of h 
would require dealing with a very large number (of the 
order of 10 6 ) of variables {pi}- This appears to be com- 
putationally intractable. 

Our phase diagram is a mean-field one - possible ef- 
fects of fluctuations are not included in our calculation. 
The first-order crystallization transition is not expected 
to be strongly affected by fluctuations. The situation is 
more complex for the glass transition because there are 
a large number of glassy local minima. When the ef- 
fects of fluctuations are included, the system might visit 
a large number of different glassy minima during its evo- 
lution over a long period of time, and thus behave like 
a liquid in the sense that the particles would no longer 
be localized in space and the time-averaged local density 
would be only weakly inhomogeneous. A true thermo- 
dynamic glass transition would occur only if the charac- 
teristic time scale for transitions between different glassy 
minima diverges in the thermodynamic limit. Whether 
this happens in the pure system is still a highly contro- 
versial issue. Further investigations of this question for 
systems of particles in the presence of quenched disorder 
would be very worthwhile. Also, the presence of multiple 
low-lying glassy minima of the free energy is expected to 
lead to slow relaxation even if no thermodynamic glass 
transition is present. Therefore, signatures of the mean- 
field glass transition found in our study should show up 
in the dynamics of the system even if no thermodynamic 
glass transition occurs when fluctuations are taken into 
consideration. 

Our density-disorder phase diagram exhibits qualita- 
tive similarities with the field-temperature phase dia- 
gram of some high-T c superconductors in the presence 
of random point pinning. For a system of vortices in the 
mixed phase of type-II superconductors, the temperature 
T plays the role of the density n* of the hard sphere sys- 
tem - increasing T is analogous to decreasing n*. As 
pointed out in the Introduction, increasing the magnetic 
field H is believed Jll[] to increase the effective strength 
of the pinning disorder. Using these analogies, one can 
translate, in a very crude and qualitative sense, our phase 
diagram in the (n* , s) plane to a phase diagram for su- 
perconductors in the (T, H) plane. Then, our result that 
the crystallization transition at weak disorder is replaced 
by a continuous glass transition as the strength of the 
disorder is increased translates into the statement that 
for superconductors, the first-order liquid to Bragg glass 
transition at low fields should change over to a contin- 
uous glass transition as the field is increased. As noted 
in the Introduction, this is precisely the behavior found 
in experiments on a family of high-T c superconductors. 
The experimentally obtained phase diagram of these su- 
perconductors also exhibits a Bragg glass to amorphous 
solid transition as the field H is increased at low temper- 
atures. This is analogous to the crystal to glass transi- 
tion found in our phase diagram as the strength of the 



12 



disorder is increased at constant density. Further evi- 
dence in support of this analogy is provided by a recent 
numerical study ]34| that suggests that the high-field, 
low-temperature phase of high-T c superconductors (the 
so-called vortex glass phase) is very similar to a struc- 
tural glass. In view of these similarities, an extension of 
our calculation to a system of pancake vortices in lay- 
ered superconductors with random point pinning, using 
the appropriate form of the free energy, would be of ob- 
vious interest. 

We are not aware of any experimentally studied sys- 
tem that provides a direct and precise physical realization 
of the model studied here. Colloidal suspensions in the 
presence of a time- independent, spatially random exter- 
nal potential (produced, for example, by suitably con- 
figured laser fields |53|) would probably provide a close 
approximation to our model. Since simple liquids with 
short range pair potentials which are strongly repulsive 
at short distances behave in many ways like a hard sphere 
liquid, our calculation is expected to apply, at least qual- 
itatively, to such systems also. 
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